Load required packages for entire analysis
dlnm.pkg <- c("tidyverse", "lubridate","xts","ggthemes","PerformanceAnalytics","reshape2","rugarch","timetk","parallel","timeSeries","tseries","data.table","ggplot2","dlnm","broom","caret","gridExtra","splines","splines2","pspline","cowplot","mgcv","spi","chron","gridGraphics","grid","pscl","MASS", "AER", "Hmisc", "MuMIn", "VGAM", "forecast", "seasonal", "plotly", "ggmap", "rgeos", "tmap", "maptools", "maps", "ggfortify", "htmltools","webshot","knitr","flexdashboard", "imager", "httr", "curl", "here")
lapply(dlnm.pkg, library, character.only=TRUE)
Load Malawi map dataset
dlnmtmp <- tempfile()
download.file("https://raw.githubusercontent.com/deusthindwa/dlnm.typhoid.nts.climate.blantyre.malawi/master/data/malawi_map.zip", destfile=dlnmtmp)
unzip(dlnmtmp, exdir = ".")
malawi.map <- rgdal::readOGR(".","malawi_map")
Subset Blantyre map
blantyre1.map <- malawi.map@data$OBJECTID >289 & malawi.map@data$OBJECTID <297
blantyre2.map <- malawi.map@data$OBJECTID >308 & malawi.map@data$OBJECTID <311
blantyre3.map <- malawi.map@data$OBJECTID >342
Convert shape file to dataframe
blantyre.map <- rbind(fortify(malawi.map[blantyre1.map,]), fortify(malawi.map[blantyre2.map,]), fortify(malawi.map[blantyre3.map,]))
blantyre.map$id <- as.integer(blantyre.map$id)
Load demographics and features dataset
blantyre.demog <- read.csv(curl("https://raw.githubusercontent.com/deusthindwa/dlnm.typhoid.nts.climate.blantyre.malawi/master/data/blantyre_demog.csv"))
map.features <- read.csv(curl("https://raw.githubusercontent.com/deusthindwa/dlnm.typhoid.nts.climate.blantyre.malawi/master/data/blantyre_features.csv"))
blantyre.demog$id <- as.integer(blantyre.demog$id)
map.all <- merge(x=blantyre.map, y=blantyre.demog, by="id", x.all=TRUE)
rm(list = ls()[grep("^blantyre", ls())])
Plot the map of Blantyre
ggplot() +
geom_polygon(data=map.all,aes(x=long,y=lat,group=group,fill=popc),colour="gray50") +
theme_classic() +
theme(axis.text.x=element_text(face="bold", size=10, color="black"), axis.text.y=element_text(face="bold", size=10, color="black")) +
labs(fill="(1998 - 2008) Population censuses") +
xlab("Longitude") +
ylab("Latitude") +
geom_point(data=map.features, aes(x=long, y=lat, shape=Geolocation, size=Geolocation), color="black") +
scale_shape_manual(values=c(17, 16, 3)) +
scale_size_manual(values=c(2,4,3)) +
theme(legend.key.height=unit(0.8, "line")) +
theme(legend.key.width=unit(0.8, "line"))

Load case dataset and subset by serovar
case <-read.csv(curl("https://raw.githubusercontent.com/deusthindwa/dlnm.typhoid.nts.climate.blantyre.malawi/master/data/case.csv"))
case$case_date <- dmy(case$case_date)
case$case_count <- c(1)
case.typhi <- subset(case, organism=="typhi")
case.iNTS <- subset(case, organism=="iNTS")
Assign 0 when no case was diagnosed
case.typhi <-aggregate(case.typhi$case_count, by=list(case.typhi$case_date), FUN=sum, na.rm=TRUE)
names(case.typhi) <- c("date", "case_count")
case.typhi <- merge(case.typhi, data.table(date=seq.Date(min(case$case_date), max(case$case_date), by="day")), by="date", all=TRUE)
case.typhi[is.na(case.typhi)]<-0
case.iNTS <-aggregate(case.iNTS$case_count, by=list(case.iNTS$case_date), FUN=sum, na.rm=TRUE)
names(case.iNTS) <- c("date", "case_count")
case.iNTS <- merge(case.iNTS, data.table(date=seq.Date(min(case$case_date), max(case$case_date), by="day")), by="date", all=TRUE)
case.iNTS[is.na(case.iNTS)]<-0
Convert to XTS and aggregate monthly
case.typhi = as.xts(case.typhi[,-1,drop=FALSE], order.by=as.Date(case.typhi[,1]))
case.iNTS = as.xts(case.iNTS[,-1,drop=FALSE], order.by=as.Date(case.iNTS[,1]))
case.typhi <- apply.monthly(case.typhi, FUN=sum)
case.iNTS <- apply.monthly(case.iNTS, FUN=sum)
Load climate dataset
climate <- read.csv(curl("https://raw.githubusercontent.com/deusthindwa/dlnm.typhoid.nts.climate.blantyre.malawi/master/data/climate.csv"))
Define rain/temp between chileka & chichiri stations
climate$climate_date <- dmy(climate$date)
climate$rainfall <- (climate$chil_r+climate$chic_r)/2
climate$temperature <- ((climate$chil_mint+climate$chil_maxt)/2+(climate$chic_mint+climate$chic_maxt)/2)/2
climate <- subset(climate, select=c(climate_date, rainfall, temperature))
climate.rain <- subset(climate, select=c(climate_date, rainfall))
climate.temp <- subset(climate, select=c(climate_date, temperature))
Convert to XTS and aggregate monthly
climate.rain <- as.xts(climate.rain[,-1, drop=FALSE], order.by=as.Date(climate.rain[,1]))
climate.temp <- as.xts(climate.temp[,-1, drop=FALSE], order.by=as.Date(climate.temp[,1]))
climate.rain <- apply.monthly(climate.rain, FUN=mean)
climate.temp <- apply.monthly(climate.temp, FUN=mean)
Create tibbles for time series decomposition
case.typhi <- tk_tbl(case.typhi, preserve_index=TRUE, rename_index="date")
case.iNTS <- tk_tbl(case.iNTS, preserve_index=TRUE, rename_index="date")
climate.rain <- tk_tbl(climate.rain, preserve_index=TRUE, rename_index="date")
climate.temp <- tk_tbl(climate.temp, preserve_index=TRUE, rename_index="date")
Seasonal-adjusted NTS cases
case.iNTS.ts <- ts(na.omit(case.iNTS$case_count), frequency=12)
trend_n <- tk_tbl(exp(seasadj(mstl(log(case.iNTS.ts)))), preserve_index=FALSE)
case.iNTS$case_count_sea <- trend_n$value
case.iNTS$case_count_sea[case.iNTS$case_count_sea < 0] <- 0
setnames(case.iNTS, old="case_count", new="case_count_obs")
trend_n <-tk_tbl(exp(mstl(log(case.iNTS.ts))))
case.iNTS$case_count_tre <- trend_n$Trend
Seasonal-adjusted typhoid cases
case.typhi.ts <- ts(na.omit(case.typhi$case_count), frequency=12)
trend_n <- tk_tbl(exp(seasadj(mstl(log(case.typhi.ts+1)))), preserve_index=FALSE)
case.typhi$case_count_sea <-trend_n$value
case.typhi$case_count_sea[case.typhi$case_count_sea < 0] <- 0
setnames(case.typhi, old="case_count", new="case_count_obs")
trend_n <- tk_tbl(exp(mstl(log(case.typhi.ts+1))))
case.typhi$case_count_tre <- trend_n$Trend
Seasonal-adjusted rain
climate.rain.ts <- ts(na.omit(climate.rain$rainfall), frequency=12)
trend_n <- tk_tbl(seasadj(mstl(climate.rain.ts)), preserve_index=FALSE)
climate.rain$rainfall_sea <- trend_n$value
climate.rain$rainfall_sea[climate.rain$rainfall_sea < 0] <- 0
setnames(climate.rain, old="rainfall", new="rainfall_obs")
trend_n <-tk_tbl(mstl(climate.rain.ts))
climate.rain$rain_tre <- trend_n$Trend
Seasonal-adjusted temperature
climate.temp.ts <- ts(na.omit(climate.temp$temperature), frequency=12)
trend_n <-tk_tbl(seasadj(mstl(climate.temp.ts)), preserve_index=FALSE)
climate.temp$temperature_sea <- trend_n$value
climate.temp$temperature_sea[climate.temp$temperature_sea < 0] <- 0
setnames(climate.temp, old="temperature", new="temperature_obs")
trend_n <- tk_tbl(mstl(climate.temp.ts))
climate.temp$temp_tre <- trend_n$Trend
rm(list=ls()[grep("^trend_n", ls())])
Plots of decomposed series
x<-case.iNTS.ts %>% mstl() %>% ggfortify:::autoplot.ts(main="A",xlab="Years (2000-2015)",size=1,colour="orange2",is.date=FALSE)+theme_bw()
y<-case.typhi.ts %>% mstl() %>% ggfortify:::autoplot.ts(main="B",xlab="Years (2000-2015)",size=1,colour="red2",is.date=FALSE)+theme_bw()
z<-climate.rain.ts %>% mstl() %>% ggfortify:::autoplot.ts(main="C",xlab="Years (2000-2015)",size=1,colour="blue2",is.date=FALSE)+theme_bw()
v<-climate.temp.ts %>% mstl() %>% ggfortify:::autoplot.ts(main="D",xlab="Years (2000-2015)",size=1,colour="green2",is.date=FALSE)+theme_bw()
grid.arrange(grobs=list(x,y,z,v), ncol=4, nrow=1)

Plots of seasonal-adjusted series
Th <- theme(plot.title=element_text(hjust=0)) +
theme(axis.title.x=element_text(size=10)) +
theme(axis.title.y=element_text(size=10)) +
theme(axis.text.x=element_text(face="bold", size=10), axis.text.y=element_text(face="bold", size=10)) +
theme(legend.key.height=unit(1, "line")) +
theme(legend.key.width=unit(1, "line"))
S1 <- ggplot(as.data.frame(case.iNTS)) +
geom_line(aes(date, case_count_obs, color="Original"), size=0.8) +
geom_line(aes(date, case_count_sea, color="Seasonal-adjusted"), size=0.8) +
scale_color_manual(values = c("Original"="black","Seasonal-adjusted"="orange2")) +
labs(title="A", x="", y="NTS cases") + theme(plot.title=element_text(hjust=0)) +
theme(legend.justification=c(0.5,0), legend.position=c(0.7,0.7), legend.text=element_text(size=10),
legend.title=element_text(size=0))
S2<-ggplot(as.data.frame(case.typhi)) +
geom_line(aes(date, case_count_obs, color="Original"), size=0.8) +
geom_line(aes(date, case_count_sea, color="Seasonal-adjusted"), size=0.8) +
scale_color_manual(values = c("Original"="black","Seasonal-adjusted"="red2")) +
labs(title="B", x ="", y="Typhoid cases") + Th +
theme(legend.justification=c(0.5,0), legend.position = c(0.3,0.7), legend.text = element_text(size = 10),
legend.title = element_text(size=0))
S3<-ggplot(as.data.frame(climate.rain)) +
geom_line(aes(date, rainfall_obs, color="Original"), size=0.8) +
geom_line(aes(date, rainfall_sea, color="Seasonal-adjusted"), size=0.8) +
scale_color_manual(values = c("Original"="black","Seasonal-adjusted"="blue2")) +
labs(title="C", x ="Year", y = "Rainfall (mm)") + Th +
theme(legend.justification=c(0.5,0), legend.position = c(0.3,0.7), legend.text = element_text(size = 10),
legend.title = element_text(size=0))
S4<-ggplot(as.data.frame(climate.temp)) +
geom_line(aes(date, temperature_obs, color="Original"), size=0.8) +
geom_line(aes(date, temperature_sea, color="Seasonal-adjusted"), size=0.8) +
scale_color_manual(values = c("Original"="black","Seasonal-adjusted"="green2")) +
labs(title="D", x ="Year", y = "Temperature (°C)") + Th +
theme(legend.justification=c(0.5,0), legend.position = c(0.3,0.7), legend.text = element_text(size = 10),
legend.title = element_text(size=0))
grid.arrange(grobs=list(S1, S2, S3, S4), ncol=2, nrow=2)

Calculate NTS/typhoid incidence
case.iNTS$census[year(case.iNTS$date)==2000]<-852054; case.iNTS$census[year(case.iNTS$date)==2001]<-873382
case.iNTS$census[year(case.iNTS$date)==2002]<-894710; case.iNTS$census[year(case.iNTS$date)==2003]<-916039
case.iNTS$census[year(case.iNTS$date)==2004]<-937367; case.iNTS$census[year(case.iNTS$date)==2005]<-958695
case.iNTS$census[year(case.iNTS$date)==2006]<-980023; case.iNTS$census[year(case.iNTS$date)==2007]<-1001352
case.iNTS$census[year(case.iNTS$date)==2008]<-1022680; case.iNTS$census[year(case.iNTS$date)==2009]<-1044008
case.iNTS$census[year(case.iNTS$date)==2010]<-1065337; case.iNTS$census[year(case.iNTS$date)==2011]<-1086665
case.iNTS$census[year(case.iNTS$date)==2012]<-1107993; case.iNTS$census[year(case.iNTS$date)==2013]<-1129322
case.iNTS$census[year(case.iNTS$date)==2014]<-1150650; case.iNTS$census[year(case.iNTS$date)==2015]<-1171978
case.iNTS$incid_sea <-case.iNTS$case_count_sea*100000/case.iNTS$census
case.iNTS$incid_obs <-case.iNTS$case_count_obs*100000/case.iNTS$census
case.typhi$census[year(case.typhi$date)==2000]<-852054; case.typhi$census[year(case.typhi$date)==2001]<-873382
case.typhi$census[year(case.typhi$date)==2002]<-894710; case.typhi$census[year(case.typhi$date)==2003]<-916039
case.typhi$census[year(case.typhi$date)==2004]<-937367; case.typhi$census[year(case.typhi$date)==2005]<-958695
case.typhi$census[year(case.typhi$date)==2006]<-980023; case.typhi$census[year(case.typhi$date)==2007]<-1001352
case.typhi$census[year(case.typhi$date)==2008]<-1022680; case.typhi$census[year(case.typhi$date)==2009]<-1044008
case.typhi$census[year(case.typhi$date)==2010]<-1065337; case.typhi$census[year(case.typhi$date)==2011]<-1086665
case.typhi$census[year(case.typhi$date)==2012]<-1107993; case.typhi$census[year(case.typhi$date)==2013]<-1129322
case.typhi$census[year(case.typhi$date)==2014]<-1150650; case.typhi$census[year(case.typhi$date)==2015]<-1171978
case.typhi$incid_sea <-case.typhi$case_count_sea*100000/case.typhi$census
case.typhi$incid_obs <-case.typhi$case_count_obs*100000/case.typhi$census
Plots of trends in NTS, typhoid, rain and temperature
E1<-ggplot() +
geom_line(data=case.iNTS, aes(x=date, y=case_count_tre, color="Number of NTS cases"), stat="identity", size=1.0) +
geom_line(data=climate.rain, aes(x=date, y=rain_tre/0.04, color="Rainfall"), alpha=0.8, size=0.7) +
scale_y_continuous(sec.axis = sec_axis(~.*0.04, name="(mm)")) + theme_bw() +
scale_color_manual(values = c("Number of NTS cases"="orange2", "Rainfall"="blue2")) +
ggtitle("A") + ylab("Cases") + xlab("Month'Year") +
theme(axis.title.x = element_text(size=0), axis.title.y=element_text(size=10, color="orange2",face="bold"),
plot.title = element_text(hjust=0,face="bold")) +
theme(axis.text.x=element_text(size=0), axis.title.y.right=element_text(color="blue2", size=10),
axis.text.y = element_text(face="bold", size=10)) +
scale_x_date(date_breaks = "24 month", date_labels ="%b'%y") +
theme(legend.justification=c(0.5,0), legend.position=c(0.6,0.7), legend.text=element_text(size=10),
legend.title=element_text(size=0)) + theme(legend.key.height=unit(1,"line")) + theme(legend.key.width=unit(1,"line"))
E2<-ggplot() +
geom_line(data=case.iNTS, aes(x=date, y=case_count_tre, color="Number of NTS cases"), stat="identity", size=1.0) +
geom_line(data=climate.temp, aes(x=date, y=temp_tre/0.4, color="Temperature"), alpha=0.8, size=0.7) +
scale_y_continuous(sec.axis = sec_axis(~.*0.4, name="(°C)")) + theme_bw() +
scale_color_manual(values=c("Number of NTS cases"="orange2", "Temperature"="green2")) +
ggtitle("B") + ylab("Cases") + xlab("Month'Year") +
theme(axis.title.x = element_text(size=10), axis.title.y=element_text(size=10, color="orange2",face="bold"),
plot.title = element_text(hjust=0,face="bold")) +
theme(axis.text.x=element_text(size=8), axis.title.y.right=element_text(color="green2", size=10),
axis.text.y=element_text(face="bold", size=10)) +
scale_x_date(date_breaks = "24 month", date_labels ="%b'%y") +
theme(legend.justification=c(0.5,0), legend.position=c(0.6,0.6), legend.text = element_text(size=10),
legend.title=element_text(size=0)) + theme(legend.key.height=unit(1,"line")) + theme(legend.key.width=unit(1,"line"))
E3<-ggplot() +
geom_line(data=case.typhi, aes(x=date, y=case_count_tre, color="Number of typhoid cases"), stat="identity", size=1.0) +
geom_line(data=climate.rain, aes(x=date, y=rain_tre/0.06, color="Rainfall"), alpha=0.8, size=0.7) +
scale_y_continuous(sec.axis=sec_axis(~.*0.06, name="(mm)")) + theme_bw() +
scale_color_manual(values=c("Number of typhoid cases"="red2", "Rainfall"="blue2")) +
ggtitle("C") + ylab("Cases") + xlab("Month'Year") +
theme(axis.title.x = element_text(size=0), axis.title.y=element_text(size=10, color="red2",face="bold"),
plot.title=element_text(hjust=0,face="bold")) +
theme(axis.text.x=element_text(size=0), axis.title.y.right=element_text(color="blue2", size=10),
axis.text.y = element_text(face="bold", size=10)) +
scale_x_date(date_breaks = "24 month", date_labels ="%b'%y") +
theme(legend.justification=c(0.5,0), legend.position=c(0.35,0.09), legend.text=element_text(size=10),
legend.title = element_text(size=0)) + theme(legend.key.height=unit(1,"line")) + theme(legend.key.width=unit(1,"line"))
E4<-ggplot() +
geom_line(data=case.typhi, aes(x=date, y=case_count_tre, color="Number of typhoid cases"), stat="identity", size=1.0) +
geom_line(data=climate.temp, aes(x=date, y=temp_tre/0.6, color="Temperature"), alpha=0.8, size = 0.7) +
scale_y_continuous(sec.axis=sec_axis(~.*0.6, name="(°C)")) + theme_bw() +
scale_color_manual(values=c("Number of typhoid cases"="red2", "Temperature"="green2")) +
ggtitle("D") + ylab("Cases") + xlab("Month'Year") +
theme(axis.title.x=element_text(size=10), axis.title.y=element_text(size=10, color="red2",face="bold"),
plot.title=element_text(hjust=0,face="bold")) +
theme(axis.text.x=element_text(face="bold", size=8), axis.title.y.right = element_text(color="green2", size=10),
axis.text.y = element_text(face="bold", size=10)) +
scale_x_date(date_breaks="24 month", date_labels="%b'%y") +
theme(legend.justification=c(0.5,0), legend.position=c(0.35,0.09), legend.text=element_text(size=10),
legend.title = element_text(size=0)) + theme(legend.key.height=unit(1,"line")) + theme(legend.key.width=unit(1,"line"))
grid.arrange(grobs=list(E1, E3, E2, E4), ncol=2, nrow=2)

Plots of Age-sex distribution of cases
case$sex[case$sex == ""] <- NA
case$age[case$age == ""] <- NA
case$date <- ymd(case$case_date)
case$year <- year(case$date)
agesex.p1 <- ggplot(subset(case, organism=="iNTS" & !is.na(age) & sex != "Unknown"), aes(x=age, color=sex)) +
geom_freqpoly(position=position_dodge(width=1.5), binwidth=1, size=1) +
scale_color_manual(values=c(Male="gray40",Female="red3")) +
theme_bw() +
scale_x_continuous(limits=c(0, 91), breaks=seq(0, 91, 5)) +
labs(title="A", x="Age (years)", y="Number of NTS cases") +
theme(plot.title=element_text(hjust=0,face="bold")) +
theme(axis.text.x=element_text(face="bold", size=10, color="black"), axis.text.y=element_text(face="bold", size=10, color="black")) +
theme(legend.justification=c(0.5,0), legend.position=c(0.5, 0.6), legend.text=element_text(size=10), legend.title=element_text(size=10)) +
labs(color="Missing sex: 2,596 (32.3%)") +
theme(legend.key.height=unit(1,"line")) +
theme(legend.key.width=unit(1,"line"))
agesex.p2<- ggplot(subset(case, organism=="typhi" & !is.na(age) & sex !="Unknown"), aes(x=age, color=sex)) +
geom_freqpoly(position=position_dodge(width=1.5), binwidth=1, size=1) +
scale_color_manual(values=c(Male="gray40",Female="red3")) +
theme_bw() +
scale_x_continuous(limits=c(0, 91), breaks=seq(0, 91, 5)) +
labs(title="B", x="Age (years)", y="Number of typhoid cases") +
theme(plot.title = element_text(hjust=0,face="bold")) +
theme(axis.text.x=element_text(face="bold", size=10, color="black"), axis.text.y=element_text(face="bold", size=10, color="black")) +
theme(legend.justification=c(0.5,0), legend.position=c(0.5, 0.6), legend.text=element_text(size=10), legend.title=element_text(size=10)) +
labs(color="Missing sex: 61 (2.4%)") +
theme(legend.key.height=unit(1,"line")) +
theme(legend.key.width=unit(1,"line"))
grid.arrange(grobs=list(agesex.p1, agesex.p2), ncol=2, nrow=1)

Generating yearly and monthly dynamics
case.iNTS.spi <- subset(case.iNTS, year(case.iNTS$date)<2011, select=c(date,incid_sea))
case.iNTS.spi$month <-month(case.iNTS.spi$date)
case.iNTS.spi$year <- year(case.iNTS.spi$date); case.iNTS.spi$date <- NULL
case.iNTS.spi <- spread(case.iNTS.spi, year, incid_sea); case.iNTS.spi <- as.matrix(case.iNTS.spi)[,-1]
YMD1 <- plot_ly(x=c(2000,2001,2002,2003,2004,2005,2006,2007,2008,2009,2010), y=c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"),
z=~case.iNTS.spi, type="contour", colorscale='jet', contours=list(showlabels=TRUE)) %>%
colorbar(title="NTS incidence per \n 100,000 population") %>%
layout(title="A", xaxis=list(title="Year"), yaxis=list(title="Month"), font=list(size=13))
case.typhi.spi <- subset(case.typhi, year(case.typhi$date)>2010, select=c(date,incid_sea))
case.typhi.spi$month <- month(case.typhi.spi$date)
case.typhi.spi$year <- year(case.typhi.spi$date); case.typhi.spi$date <- NULL
case.typhi.spi <- spread(case.typhi.spi, year, incid_sea); case.typhi.spi <- as.matrix(case.typhi.spi)[,-1]
YMD2 <- plot_ly(x = c(2011,2012,2013,2014,2015), y=c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"),
z=~case.typhi.spi, type="contour", colorscale='heatmap', contours=list(showlabels=TRUE)) %>%
colorbar(title="Typhoid incidence per \n 100,000 population") %>%
layout(title="D", xaxis=list(title="Year"), yaxis=list(title="Month"), font=list(size=13))
climate.rain.spin <- subset(climate.rain, year(climate.rain$date)<2011, select=c(date,rainfall_obs))
climate.rain.spin$month <- month(climate.rain.spin$date)
climate.rain.spin$year <- year(climate.rain.spin$date); climate.rain.spin$date <- NULL
climate.rain.spin <- spread(climate.rain.spin, year, rainfall_obs); climate.rain.spin <- as.matrix(climate.rain.spin)[,-1]
YMD3 <- plot_ly(x = c(2000,2001,2002,2003,2004,2005,2006,2007,2008,2009,2010), y=c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"),
z=~climate.rain.spin, type="contour", colorscale='heatmap', contours=list(showlabels=TRUE)) %>% colorbar(title="Rainfall (mm)") %>%
layout(title="B", xaxis=list(title="Year"), yaxis=list(title="Month"), font=list(size=13))
climate.rain.spit <- subset(climate.rain, year(climate.rain$date)>2010, select=c(date,rainfall_obs))
climate.rain.spit$month <- month(climate.rain.spit$date)
climate.rain.spit$year <- year(climate.rain.spit$date); climate.rain.spit$date <- NULL
climate.rain.spit <- spread(climate.rain.spit, year, rainfall_obs); climate.rain.spit <- as.matrix(climate.rain.spit)[,-1]
YMD4 <- plot_ly(x = c(2011,2012,2013,2014,2015), y=c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"),
z=~climate.rain.spit, type="contour", colorscale='heatmap', contours=list(showlabels=TRUE)) %>%
colorbar(title="Rainfall (mm)") %>%
layout(title="E", xaxis=list(title="Year"), yaxis=list(title="Month"), font=list(size=13))
climate.temp.spin <- subset(climate.temp, year(climate.temp$date)<2011, select=c(date,temperature_obs))
climate.temp.spin$month <- month(climate.temp.spin$date)
climate.temp.spin$year <- year(climate.temp.spin$date); climate.temp.spin$date <- NULL
climate.temp.spin <- spread(climate.temp.spin, year, temperature_obs); climate.temp.spin <- as.matrix(climate.temp.spin)[,-1]
YMD5 <- plot_ly(x=c(2000,2001,2002,2003,2004,2005,2006,2007,2008,2009,2010), y=c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"),
z=~climate.temp.spin, type="contour", colorscale='heatmap', contours=list(showlabels=TRUE)) %>%
colorbar(title="Temperature (°C)") %>%
layout(title="C", xaxis=list(title="Year"), yaxis=list(title="Month"), font=list(size=13))
climate.temp.spit <- subset(climate.temp, year(climate.temp$date)>2010, select=c(date,temperature_obs))
climate.temp.spit$month <- month(climate.temp.spit$date)
climate.temp.spit$year <- year(climate.temp.spit$date); climate.temp.spit$date <- NULL
climate.temp.spit <- spread(climate.temp.spit, year, temperature_obs); climate.temp.spit <- as.matrix(climate.temp.spit)[,-1]
YMD6 <- plot_ly(x=c(2011,2012,2013,2014,2015), y=c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"),
z=~climate.temp.spit, type="contour", colorscale='heatmap', contours=list(showlabels=TRUE)) %>%
colorbar(title="Temperature (°C)") %>%
layout(title="F", xaxis=list(title="Year"), yaxis=list(title="Month"), font=list(size = 13))
Plots of yearly and monthly dynamics
htmltools::tags$div(
style = "display: flex; flex-wrap: wrap",
tags$div(YMD1, style = "width: 33%; padding: 1em;"),
tags$div(YMD3, style = "width: 33%; padding: 1em;"),
tags$div(YMD5, style = "width: 33%; padding: 1em;"),
tags$div(YMD2, style = "width: 33%; padding: 1em;"),
tags$div(YMD4, style = "width: 33%; padding: 1em;"),
tags$div(YMD6, style = "width: 33%; padding: 1em;")
)
Constructing DLNM datasets for NTS and Typhoid
mo.dlnmN <- bind_cols(case.iNTS, climate.rain, climate.temp, id=NULL)
mo.dlnmN$date1 <- mo.dlnmN$date2 <- NULL
mo.dlnmN <- subset(mo.dlnmN, year(date) < 2011)
mo.dlnmN$time <- seq.int(from = 1, to=132, by=1)
mo.dlnmN$year <- year(mo.dlnmN$date)
mo.dlnmN$month <- month(mo.dlnmN$date)
mo.dlnmN$incid_seaX<-round(mo.dlnmN$incid_sea, digits = 0)
mo.dlnmN$incid_obsX<-round(mo.dlnmN$incid_obs, digits = 0)
mo.dlnmT <- bind_cols(case.typhi, climate.rain, climate.temp, id=NULL)
mo.dlnmT$date1 <- mo.dlnmT$date2 <- NULL
mo.dlnmT <- subset(mo.dlnmT, year(date) > 2010)
mo.dlnmT$time <- seq.int(from = 1, to=60, by=1)
mo.dlnmT$year <- year(mo.dlnmT$date)
mo.dlnmT$month <- month(mo.dlnmT$date)
mo.dlnmT$incid_seaX<-round(mo.dlnmT$incid_sea, digits = 0)
AIC to QAIC for NTS and Typhoid
nts.quasipoisson <- function(...) {
res <- quasipoisson(...)
res$aic <- poisson(...)$aic
res
}
typ.quasipoisson <- function(...) {
res <- quasipoisson(...)
res$aic <- poisson(...)$aic
res
}
Calculate QAIC values to select optimal models
QAICtable <- data.frame(model.no=rep(NA,27), lag.df=rep(NA,27), fx1.df=rep(NA,27), fx2.df=rep(NA,27),
QAIC.ntsR=rep(NA,27), QAIC.ntsT=rep(NA,27), QAIC.nts=rep(NA,27), QAIC.typR=rep(NA,27),
QAIC.typT=rep(NA,27), QAIC.typ=rep(NA,27))
l=1
for(k in 3:5){
for(j in 3:5){
for(i in 3:5){
nts.lagknots <- logknots(8, fun="ns", df=i)
nts.varknots.r=equalknots(mo.dlnmN$rainfall_obs, fun="ns", df=j)
nts.varknots.t=equalknots(mo.dlnmN$temperature_obs, fun="ns", df=k)
nts.mo.cb.rain <- crossbasis(mo.dlnmN$rainfall_obs, lag=8, argvar=list(fun="ns", knots=nts.varknots.r), arglag=list(knots=nts.lagknots))
nts.mo.cb.temp <- crossbasis(mo.dlnmN$temperature_obs, lag=8, argvar=list(fun="ns", knots=nts.varknots.t), arglag=list(knots=nts.lagknots))
nts.modelR <- glm(mo.dlnmN$incid_seaX ~ nts.mo.cb.rain + year, family=nts.quasipoisson(), na.action=na.delete, mo.dlnmN)
nts.modelT <- glm(mo.dlnmN$incid_seaX ~ nts.mo.cb.temp + year, family=nts.quasipoisson(), na.action=na.delete, mo.dlnmN)
nts.model <- glm(mo.dlnmN$incid_seaX ~ nts.mo.cb.rain + nts.mo.cb.temp + year, family = nts.quasipoisson(), na.action=na.delete, mo.dlnmN)
typ.lagknots <- logknots(8, fun="ns", df=i)
typ.varknots.r=equalknots(mo.dlnmT$rainfall_obs, fun="ns", df=j)
typ.varknots.t=equalknots(mo.dlnmT$temperature_obs, fun="ns", df=k)
typ.mo.cb.rain <- crossbasis(mo.dlnmT$rainfall_obs, lag=8, argvar=list(fun="ns", knots=typ.varknots.r), arglag=list(knots=typ.lagknots))
typ.mo.cb.temp <- crossbasis(mo.dlnmT$temperature_obs, lag =8, argvar = list(fun="ns", knots=typ.varknots.t), arglag=list(knots=typ.lagknots))
typ.modelR <- glm(mo.dlnmT$incid_seaX ~ typ.mo.cb.rain + year, family = typ.quasipoisson(), na.action=na.delete, mo.dlnmT)
typ.modelT <- glm(mo.dlnmT$incid_seaX ~ typ.mo.cb.temp + year, family = typ.quasipoisson(), na.action=na.delete, mo.dlnmT)
typ.model <- glm(mo.dlnmT$incid_seaX ~ typ.mo.cb.rain + typ.mo.cb.temp + year, family = typ.quasipoisson(), na.action=na.delete, mo.dlnmT)
QAICtable[l,] <- c(l,i,j,k, QAIC(nts.modelR, chat=summary(nts.modelR)$dispersion), QAIC(nts.modelT, chat=summary(nts.modelT)$dispersion),
QAIC(nts.model, chat=summary(nts.model)$dispersion), QAIC(typ.modelR, chat=summary(typ.modelR)$dispersion), QAIC(typ.modelT, chat=summary(typ.modelT)$dispersion),
QAIC(typ.model, chat=1))
l=l+1
}
}
}
kable(QAICtable)
| 1 |
3 |
3 |
3 |
570.8194 |
565.8631 |
562.3950 |
218.5004 |
220.2006 |
212.6658 |
| 2 |
4 |
3 |
3 |
576.3831 |
571.4838 |
573.3630 |
223.1714 |
216.1853 |
220.6390 |
| 3 |
5 |
3 |
3 |
581.4852 |
574.9872 |
583.4617 |
228.2399 |
221.1698 |
230.8468 |
| 4 |
3 |
4 |
3 |
571.1280 |
565.8631 |
565.3731 |
211.3436 |
220.2006 |
216.2219 |
| 5 |
4 |
4 |
3 |
578.3398 |
571.4838 |
578.1025 |
218.1843 |
216.1853 |
228.0142 |
| 6 |
5 |
4 |
3 |
585.4681 |
574.9872 |
590.3358 |
222.7394 |
221.1698 |
237.6362 |
| 7 |
3 |
5 |
3 |
572.7386 |
565.8631 |
566.1964 |
216.3176 |
220.2006 |
221.5570 |
| 8 |
4 |
5 |
3 |
581.6967 |
571.4838 |
581.3661 |
222.4327 |
216.1853 |
234.0871 |
| 9 |
5 |
5 |
3 |
590.5606 |
574.9872 |
594.8792 |
227.8628 |
221.1698 |
246.2023 |
| 10 |
3 |
3 |
4 |
570.8194 |
571.7243 |
566.2250 |
218.5004 |
220.7893 |
215.5219 |
| 11 |
4 |
3 |
4 |
576.3831 |
577.5401 |
578.9833 |
223.1714 |
220.2852 |
225.9024 |
| 12 |
5 |
3 |
4 |
581.4852 |
583.4447 |
589.9669 |
228.2399 |
226.2618 |
236.0380 |
| 13 |
3 |
4 |
4 |
571.1280 |
571.7243 |
568.6740 |
211.3436 |
220.7893 |
220.0702 |
| 14 |
4 |
4 |
4 |
578.3398 |
577.5401 |
583.3435 |
218.1843 |
220.2852 |
233.1946 |
| 15 |
5 |
4 |
4 |
585.4681 |
583.4447 |
597.5272 |
222.7394 |
226.2618 |
245.6773 |
| 16 |
3 |
5 |
4 |
572.7386 |
571.7243 |
571.2271 |
216.3176 |
220.7893 |
224.5520 |
| 17 |
4 |
5 |
4 |
581.6967 |
577.5401 |
588.2408 |
222.4327 |
220.2852 |
239.9495 |
| 18 |
5 |
5 |
4 |
590.5606 |
583.4447 |
603.1850 |
227.8628 |
226.2618 |
252.6551 |
| 19 |
3 |
3 |
5 |
570.8194 |
574.0369 |
566.7772 |
218.5004 |
220.7423 |
220.4298 |
| 20 |
4 |
3 |
5 |
576.3831 |
580.7484 |
580.7418 |
223.1714 |
220.9285 |
232.7209 |
| 21 |
5 |
3 |
5 |
581.4852 |
587.7540 |
593.8723 |
228.2399 |
225.5709 |
245.1746 |
| 22 |
3 |
4 |
5 |
571.1280 |
574.0369 |
569.9469 |
211.3436 |
220.7423 |
224.8192 |
| 23 |
4 |
4 |
5 |
578.3398 |
580.7484 |
586.2703 |
218.1843 |
220.9285 |
239.8864 |
| 24 |
5 |
4 |
5 |
585.4681 |
587.7540 |
601.6447 |
222.7394 |
225.5709 |
254.6679 |
| 25 |
3 |
5 |
5 |
572.7386 |
574.0369 |
575.0626 |
216.3176 |
220.7423 |
229.8974 |
| 26 |
4 |
5 |
5 |
581.6967 |
580.7484 |
593.7784 |
222.4327 |
220.9285 |
246.7405 |
| 27 |
5 |
5 |
5 |
590.5606 |
587.7540 |
610.6531 |
227.8628 |
225.5709 |
261.6885 |
Fit NTS model
nts.model <- glm(mo.dlnmN$incid_seaX ~ nts.cb.rain + nts.cb.temp + year, family=quasipoisson(), na.action=na.delete, mo.dlnmN)
Investigate deviance residuals for a fitted NTS model
par(mfrow=c(1,2))
pacf(residuals(nts.model,type="deviance"), na.action=na.omit, main="Partial autocorrelation (original modelR+T)", xlim=c(0,8))
nts.model <- update(nts.model,.~.+Lag(residuals(nts.model,type="deviance"), 1)) #add residuals at lag 1
pacf(residuals(nts.model,type="deviance"), na.action=na.omit, main="Partial autocorrelation (adjusted modelR+T)", xlim=c(0,8))

Run predictions
nts.pred.rain <- crosspred(nts.cb.rain, nts.model, cen=0, by=0.2)
nts.pred.temp <- crosspred(nts.cb.temp, nts.model, cen=23, by=0.2)
Plots of NTS predictions due to rainfall
plot(nts.pred.rain, xlab="Rainfall (mm)", ylab="Lag (month)", zlab="RR of NTS", theta=150, phi=5, lphi=100, cex.lab=1.1,
cex.axis=1.1, col="gray80",main="A")

plot(nts.pred.rain, "contour", key.title=title("NTS"), plot.title=title("B", xlab ="Rainfall (mm)", ylab = "Lag (month)",
cex.lab=1.1, cex.axis=1.1))

plot(nts.pred.rain, "slices", xlab="Rainfall (mm)", lag=c(2), col="orange2", ci.arg=list(col=topo.colors(70, alpha = 1)),
ci.level=0.95, ci='b', lwd=4.5, ylab="RR of NTS", cex.lab=1.1, cex.axis=1.1,main="C")

plot(nts.pred.rain, "slices", xlab="Lag (month)", var=c(9), col="orange2", ci.arg=list(col=topo.colors(70, alpha = 1)),
ci.level=0.95, ci='b', lwd=4.5, ylab="RR of NTS", cex.lab=1.1, cex.axis=1.1,main="D")

Plots of NTS predictions due to temperature
plot(nts.pred.temp, xlab="Temperature (°C)", ylab="Lag (month)", zlab="RR of NTS", theta=150, phi=5, lphi=100, cex.lab=1, cex.axis=1, col="gray80",main="A")

plot(nts.pred.temp, "contour", key.title=title("NTS"), plot.title=title("B", xlab ="Temperature (°C)", ylab = "Lag (month)", cex.lab=1, cex.axis=1))

plot(nts.pred.temp, xlab="Temperature (°C)", "slices",lag=c(3), col="orange2", ci.arg=list(col=terrain.colors(70, alpha = 1)), ci.level=0.95, ci='b',lwd=4.5, ylab="RR of NTS", cex.lab=1.1, cex.axis=1.1,main="C")

plot(nts.pred.temp, xlab="Lag (month)", "slices",var=c(19), col="orange2", ci.arg=list(col=terrain.colors(70, alpha = 1)), ci.level=0.95, ci='b',lwd=4.5, ylab="RR of NTS", cex.lab=1.1, cex.axis=1.1,main="D")

Fit Typhoid model
typ.modelR <- glm(mo.dlnmT$incid_seaX ~ typ.cb.rain1 + year, family=quasipoisson(), na.action=na.delete, mo.dlnmT)
typ.modelT <- glm(mo.dlnmT$incid_seaX ~ typ.cb.temp + typ.cb.rain2 + year, family=quasipoisson(), na.action=na.delete, mo.dlnmT)
Investigate deviance residuals for fitted tyhpoid models
par(mfrow=c(1,3))
pacf(residuals(typ.modelR,type="deviance"),na.action=na.omit,main="Partial Autocorrelation (original modelR)",xlim=c(0,8))
pacf(residuals(typ.modelT,type="deviance"),na.action=na.omit,main="Partial Autocorrelation (original modelT)",xlim=c(0,8))
typ.modelT <- update(typ.modelT,.~.+Lag(residuals(typ.modelT,type="deviance"),1)) #add residuals at lag 1
pacf(residuals(typ.modelT,type="deviance"),na.action=na.omit,main="Partial Autocorrelation (adjusted modelT)",xlim=c(0,8))

Run predictions
typ.pred.rain <- crosspred(typ.cb.rain1, typ.modelR, cen=0, by=0.2)
typ.pred.temp <- crosspred(typ.cb.temp, typ.modelT, cen=23, by=0.2)
Plots of TYP predictions due to rainfall
plot(typ.pred.rain, xlab="Rainfall (mm)", ylab="Lag (month)", zlab="RR of typhoid", theta=150, phi=5, lphi=150, col="gray80", main="A", cex.lab=1.1, cex.axis=1.1)

plot(typ.pred.rain, "contour", key.title=title("TYP"), plot.title=title("B", xlab ="Rainfall (mm)", ylab = "Lag (month)", cex.lab=1.1, cex.axis=1.1))

plot(typ.pred.rain, xlab="Rainfall (mm)", "slices", lag=c(4), col="red2", ci.arg=list(col=topo.colors(70, alpha = 1)), ci.level=0.95,ci='b', lwd=4.5, ylab="RR of typhoid", cex.lab=1.1, cex.axis=1.1,main="C")

plot(typ.pred.rain, xlab="Lag (month)", "slices", var=c(9), col="red2", ci.arg=list(col=topo.colors(70, alpha = 1)), ci.level=0.95,ci='b', lwd=4.5, ylab="RR of typhoid", cex.lab=1.1, cex.axis=1.1,main="D")

Plots of TYP predictions due to temperature
plot(typ.pred.temp, xlab="Temperature (°C)", ylab="Lag (month)", zlab="RR of typhoid", theta=150, phi=5, lphi=100, col="gray80", main="A", cex.lab=1.1, cex.axis=1.1)

plot(typ.pred.temp, "contour", key.title=title("TYP"), plot.title=title("B", xlab ="Temperature (°C)", ylab = "Lag (month)", cex.lab=1.1, cex.axis=1.1))

plot(typ.pred.temp, xlab="Temperature (°C)", "slices", lag=c(2), col="red2", ci.arg=list(col=terrain.colors(70, alpha = 1)), ci.level=0.95,ci='b', lwd=4.5, ylab="RR of typhoid", cex.lab=1.1, cex.axis=1.1,main="C")

plot(typ.pred.temp, xlab="Temperature (°C)", "slices", lag=c(5), col="red2", ci.arg=list(col=terrain.colors(70, alpha = 1)), ci.level=0.95,ci='b', lwd=4.5, ylab="RR of typhoid", cex.lab=1.1, cex.axis=1.1,main="D")

plot(typ.pred.temp, xlab="Lag (month)", "slices", var=c(19), col="red2", ci.arg=list(col=terrain.colors(70, alpha = 1)), ci.level=0.95,ci='b', lwd=4.5, ylab="RR of typhoid", cex.lab=1.1, cex.axis=1.1,main="E")

plot(typ.pred.temp, xlab="Lag (month)", "slices", var=c(25), col="red2", ci.arg=list(col=terrain.colors(70, alpha = 1)), ci.level=0.95,ci='b', lwd=4.5, ylab="RR of typhoid", cex.lab=1.1, cex.axis=1.1,main="F")

Plots of validation checks
par(mfrow=c(3,4))
plot(mo.dlnmN$date, residuals(nts.model,type="deviance"), pch=19, cex=0.8, col=grey(0.4),main="A", ylab="Residuals",xlab="",cex.lab=1.2,cex.axis=1.2)
abline(h=0,lty=2,lwd=3)
pacf(residuals(nts.model,type="deviance"),na.action=na.omit,main="B",xlim=c(0,8),xlab="", cex.lab=1.2,cex.axis=1.2)
acf(residuals(nts.model,type="deviance"),na.action=na.omit,main="C",xlim=c(0,8),xlab="", cex.lab=1.2,cex.axis=1.2)
plot(mo.dlnmN$incid_seaX, pch=10, cex=0.8, col=grey(0.6), size =1 ,main="D", ylab="NTS incidence",xlab="",cex.lab=1.2, cex.axis=1.2)
lines(predict(nts.model, type="response"), col="orange2",lwd=3)
legend("topright", legend=c("observed", "predicted"), col=c("grey", "orange2"), lty=1:1, cex=0.8, lwd=3)
plot(mo.dlnmT$date, residuals(typ.modelR,type="deviance"), pch=19, cex=0.8, col=grey(0.6),main="E", ylab="Residuals",xlab="",cex.lab=1.2,cex.axis=1.2)
abline(h=0,lty=2,lwd=3)
pacf(residuals(typ.modelR,type="deviance"),na.action=na.omit,main="F",xlim=c(0,8),xlab="",cex.lab=1.2, cex.axis=1.2)
acf(residuals(typ.modelR,type="deviance"),na.action=na.omit,main="G",xlim=c(0,8),xlab="",cex.lab=1.2, cex.axis=1.2)
plot(mo.dlnmT$incid_seaX, pch=10, cex=0.8, col=grey(0.6), size =1 ,main="H", ylab="Typhoid incidence",xlab="",cex.lab=1.2, cex.axis=1.2)
lines(predict(typ.modelR, type="response"), col="red",lwd=3)
legend("topleft", legend=c("observed", "predicted"), col=c("grey", "red2"), lty=1:1, cex=0.8, lwd=3)
plot(mo.dlnmT$date, residuals(typ.modelT,type="deviance"), pch=19, cex=0.8, col=grey(0.6),main="I", ylab="Residuals", xlab="Year",cex.lab=1.2,cex.axis=1.2)
abline(h=0,lty=2,lwd=3)
pacf(residuals(typ.modelT,type="deviance"),na.action=na.omit,main="J",xlim=c(0,8),xlab="Lag (month)",cex.lab=1.2, cex.axis=1.2)
acf(residuals(typ.modelT,type="deviance"),na.action=na.omit,main="K",xlim=c(0,8),xlab="Lag (month)",cex.lab=1.2, cex.axis=1.2)
plot(mo.dlnmT$incid_seaX, pch=10, cex=0.8, col=grey(0.6), size =1 ,main="L", ylab="Typhoid incidence", xlab="Month",cex.lab=1.2, cex.axis=1.2)
lines(predict(typ.modelT, type="response"), col="red2",lwd=3)
legend("topleft", legend=c("observed", "predicted"), col=c("grey", "red"), lty=1:1, cex=0.8, lwd=3)

Plots of TYP predictions due to temperature
cbind(nts.pred.rain$matRRfit, nts.pred.rain$matRRlow, nts.pred.rain$matRRhigh)["9",]
## lag0 lag1 lag2 lag3 lag4 lag5 lag6
## 1.2248467 1.2514626 1.2554780 1.2238951 1.1635420 1.0841819 0.9951394
## lag7 lag8 lag0 lag1 lag2 lag3 lag4
## 0.9042881 0.8176177 1.0254867 1.1162267 1.1453653 1.1143857 1.0596817
## lag5 lag6 lag7 lag8 lag0 lag1 lag2
## 0.9906237 0.9027750 0.7975713 0.6897005 1.4629635 1.4030831 1.3761766
## lag3 lag4 lag5 lag6 lag7 lag8
## 1.3441658 1.2775817 1.1865762 1.0969538 1.0252839 0.9692593
cbind(nts.pred.temp$matRRfit, nts.pred.temp$matRRlow, nts.pred.temp$matRRhigh)["19",]
## lag0 lag1 lag2 lag3 lag4 lag5 lag6
## 1.0822001 1.1309828 1.1629884 1.1659357 1.1432279 1.1012306 1.0467408
## lag7 lag8 lag0 lag1 lag2 lag3 lag4
## 0.9861515 0.9249535 0.9602616 1.0318836 1.0674043 1.0690932 1.0514755
## lag5 lag6 lag7 lag8 lag0 lag1 lag2
## 1.0192753 0.9708904 0.9049599 0.8287981 1.2196229 1.2395991 1.2671318
## lag3 lag4 lag5 lag6 lag7 lag8
## 1.2715505 1.2429867 1.1897756 1.1285169 1.0746274 1.0322647
cbind(nts.pred.temp$matRRfit, nts.pred.temp$matRRlow, nts.pred.temp$matRRhigh)["25",]
## lag0 lag1 lag2 lag3 lag4 lag5 lag6
## 0.9983030 0.9778710 0.9641177 0.9603166 0.9651159 0.9768977 0.9941390
## lag7 lag8 lag0 lag1 lag2 lag3 lag4
## 1.0153070 1.0387806 0.9057049 0.9118456 0.9084918 0.9062719 0.9129857
## lag5 lag6 lag7 lag8 lag0 lag1 lag2
## 0.9274517 0.9428357 0.9521165 0.9547300 1.1003682 1.0486771 1.0231495
## lag3 lag4 lag5 lag6 lag7 lag8
## 1.0175842 1.0202227 1.0289798 1.0482338 1.0826914 1.1302307
cbind(typ.pred.rain$matRRfit, typ.pred.rain$matRRlow, typ.pred.rain$matRRhigh)["9",]
## lag0 lag1 lag2 lag3 lag4 lag5 lag6
## 0.9870340 1.2979004 1.6063110 1.8078055 1.8722857 1.8143089 1.6725894
## lag7 lag8 lag0 lag1 lag2 lag3 lag4
## 1.4915122 1.3081090 0.6284634 0.9438326 1.2552773 1.4434731 1.5069006
## lag5 lag6 lag7 lag8 lag0 lag1 lag2
## 1.4634299 1.3249680 1.1234975 0.9124989 1.5501876 1.7847926 2.0555100
## lag3 lag4 lag5 lag6 lag7 lag8
## 2.2640954 2.3262675 2.2493164 2.1114137 1.9800744 1.8752341
cbind(typ.pred.temp$matRRfit, typ.pred.temp$matRRlow, typ.pred.temp$matRRhigh)["19",]
## lag0 lag1 lag2 lag3 lag4 lag5 lag6
## 1.4960476 1.5619360 1.5874626 1.5468052 1.4526163 1.3245023 1.1812571
## lag7 lag8 lag0 lag1 lag2 lag3 lag4
## 1.0380766 0.9055470 1.0424914 1.1457962 1.0968536 1.0279449 0.9761052
## lag5 lag6 lag7 lag8 lag0 lag1 lag2
## 0.9252186 0.8373032 0.6896236 0.5229565 2.1469324 2.1292129 2.2975149
## lag3 lag4 lag5 lag6 lag7 lag8
## 2.3275628 2.1617487 1.8960993 1.6665030 1.5625959 1.5680377
cbind(typ.pred.temp$matRRfit, typ.pred.temp$matRRlow, typ.pred.temp$matRRhigh)["25",]
## lag0 lag1 lag2 lag3 lag4 lag5 lag6
## 0.6881750 0.9207072 1.1640344 1.3467161 1.4416997 1.4504554 1.3928588
## lag7 lag8 lag0 lag1 lag2 lag3 lag4
## 1.2966573 1.1885051 0.4844865 0.7123388 0.8891649 0.9975957 1.0683040
## lag5 lag6 lag7 lag8 lag0 lag1 lag2
## 1.1081654 1.1070709 1.0400232 0.9069586 0.9774984 1.1900260 1.5238750
## lag3 lag4 lag5 lag6 lag7 lag8
## 1.8180153 1.9456054 1.8984720 1.7524222 1.6166181 1.5574519